Source code for BigDFT.LRTDDFT


import numpy as np
from futile.Utils import write
HaeV = 27.21138386


def _occ_and_virt(log):
    """
    Extract the number of occupied and empty orbitals from a logfile
    """
    norb = log.log['Total Number of Orbitals']
    if log.log['Spin treatment'] == 'Averaged':
        norbv = log.evals[0].info[0]-norb
        return (norb,), (norbv,)
    elif log.log['Spin treatment'] == 'Collinear':
        mpol = log.log['dft']['mpol']
        norbu = int((norb+mpol)/2)
        norbd = norb-norbu
        norbvu = log.evals[0].info[0]-norbu
        norbvd = log.evals[0].info[0]-norbd
        return (norbu, norbd), (norbvu, norbvd)
    else:
        raise ValueError('Information for the orbitals to be implemented')


[docs]def transition_indexes(np, nalpha, indexes): """ Returns the list of the indices in the bigdft convention that correspond to the couple iorb-ialpha with given spin. Args: np (tuple): (norbu,norbd) occupied orbitals: when of length 1 assumed spin averaged nalpha (tuple): (norbu, norbd)virtual orbitals: when of length 1 assumed spin averaged indexes (list): list of tuples of (iorb,ialpha,ispin) desired indices in python convention (start from 0) """ nspin = len(np) inds = [] for iorb, ialpha, ispin in indexes: jspin = ispin if nspin == 2 else 0 ind = ialpha+iorb*nalpha[jspin] # local index of the spin subspace if ispin == 1: ind += np[0]*nalpha[0] # spin 2 comes after spin one inds.append(ind) return inds
def _collection_indexes(np, nvirt_small): harvest = [] for ispin in [0, 1]: jspin = ispin if len(np) == 2 else 0 for ip in range(np[jspin]): for ialpha in range(nvirt_small[jspin]): harvest.append([ip, ialpha, ispin]) return harvest def _collection_indexes_iocc(iocc, nvirt, spin=None): """ For each iocc and a selected spin provide the list of couples that are concerned up to nvirt If spin is none provide the list for all values of the spin """ harvest = [] for ispin in [0, 1]: jspin = ispin if len(nvirt) == 2 else 0 if spin is not None and ispin != spin: continue for ialpha in range(nvirt[jspin]): harvest.append([iocc, ialpha, ispin]) return harvest
[docs]class TransitionMatrix(np.matrix): """ Matrix of Transition Quantities, might be either :class:`CouplingMatrix` or :class:`TransitionMultipoles` Args: matrix (matrix-like): data of the coupling matrix. If present also the number of orbitals should be provided. norb_occ (tuple): number of occupied orbitals per spin channnel. Compulsory if ``matrix`` is specified. norb_virt (tuple): number of empty orbitals per spin channnel. Compulsory if ``matrix`` is specified. log (Logfile): Instance of the logfile from which the coupling matrix calculation is performed. Automatically retrieves the ``norb_occ`` and `norb_virt`` parameters. When ``log`` parameter is present the parameter ``matrix`` is ignored. Raises: ValueError: if the file of the coupling matrix indicated by ``log`` does not exists """ def __new__(cls, matrix=None, norb_occ=None, norb_virt=None, log=None): """ Create the object from the arguments and return the ``self`` instance """ import os if log is not None: datadir = log.log.get('radical', '') datadir = 'data-'+datadir if len(datadir) > 0 else 'data' cmfile = os.path.join(log.srcdir, datadir, cls._filename) if not os.path.isfile(cmfile): raise ValueError('The file "'+cmfile+'" does not exist') norb, norbv = _occ_and_virt(log) write('Loading data with ', norb, ' occupied and ', norbv, ' empty states, from file "', cmfile, '"') try: import pandas as pd write('Using pandas:') mat = pd.read_csv(cmfile, delim_whitespace=True, header=None) except ImportError: write('Using numpy:') mat = np.loadtxt(cmfile) write('done') else: mat = matrix return super(TransitionMatrix, cls).__new__(cls, mat) def __init__(self, *args, **kwargs): """ Perform sanity checks on the loaded matrix """ log = kwargs.get('log') if log is not None: self.norb_occ, self.norb_virt = _occ_and_virt(log) else: self.norb_occ = kwargs.get('norb_occ') self.norb_virt = kwargs.get('norb_virt') assert(self.shape[0] == self._total_transitions()) write("Shape is conformal with the number of orbitals") self._sanity_check() def _total_transitions(self): ntot = 0 for no, nv in zip(self.norb_occ, self.norb_virt): ntot += no*nv if len(self.norb_occ) == 1: ntot *= 2 return ntot def _subindices(self, norb_occ, norb_virt): for i, (no, nv) in enumerate(zip(norb_occ, norb_virt)): assert(no <= self.norb_occ[i] and nv <= self.norb_virt[i]) harvest = _collection_indexes(norb_occ, norb_virt) return np.array(transition_indexes(norb_occ, self.norb_virt, harvest)) def _sanity_check(self): pass
[docs]class CouplingMatrix(TransitionMatrix): """ Casida Coupling Matrix, extracted from the calculation performed by BigDFT """ _filename = 'coupling_matrix.txt' def _sanity_check(self): write('Casida Matrix is symmetric', np.allclose(self, self.T, atol=1.e-12))
[docs] def subportion(self, norb_occ, norb_virt): """Extract a subportion of the coupling matrix. Returns a Coupling Matrix which is made by only considering the first ``norb_occ`` and ``norb_virt`` orbitals Args: norb_occ (tuple): new number of occupied orbitals. Must be lower that the instance value norb_virt (tuple): new number of virtual orbitals. Must be lower that the instance value """ inds = self._subindices(norb_occ, norb_virt) mat = np.array([row[0, inds] for row in self[inds]]) return CouplingMatrix(matrix=mat, norb_occ=norb_occ, norb_virt=norb_virt)
[docs] def diagonalize(self): """ Diagonalize the Coupling Matrix Returns: (np.matrix, np.matrix): tuple of the Eigenvvalues and Eigenvectors of the coupling matrix, as returned by :meth:`numpy.linalg.eigh`. We perform the transpose of the matrix with eigenvectors to have them sorted as row vectors """ write('Diagonalizing Coupling matrix of shape', self.shape) E2, C_E2 = np.linalg.eigh(self) write('Eigensystem solved') C_E2 = C_E2.T return E2, C_E2
[docs]class TransitionMultipoles(TransitionMatrix): """ Transition dipoles, extracted from the calculation performed by BigDFT """ _filename = 'transition_quantities.txt'
[docs] def subportion(self, norb_occ, norb_virt): """Extract a subportion of the Transition Multipoles. Returns a set of transition multipoles which is made by only considering the first ``norb_occ`` and ``norb_virt`` orbitals Args: norb_occ (tuple): new number of occupied orbitals. Must be lower that the instance value norb_virt (tuple): new number of virtual orbitals. Must be lower that the instance value Returns: TransitionMultipoles: reduced transition multipoles """ inds = self._subindices(norb_occ, norb_virt) mat = np.array(self[inds]) return TransitionMultipoles(matrix=mat, norb_occ=norb_occ, norb_virt=norb_virt)
[docs] def get_transitions(self): """ Get the transition quantities as the dimensional objects which should contribute to the oscillator strengths. Returns: numpy.array: Transition quantities multiplied by the square root of the unperturbed transition energy """ newdipole = [] for line in self: newdipole.append(np.ravel(line[0, 0]*line[0, 1:])) return np.array(newdipole)
[docs]class TransitionDipoles(TransitionMultipoles): """ Transition dipoles as provided in the version of the code < 1.8.0. Deprecated, to be used in some particular cases """ _filename = 'transition_dipoles.txt'
[docs] def get_transitions(self): return self
[docs]class Excitations(): """LR Excited states of a system Definition of the excited states in the Casida Formalism Args: cm (CouplingMatrix): the matrix of coupling among transitions tm (TransitionMultipoles): scalar product of multipoles among transitions """ def __init__(self, cm, tm): self.cm = cm self.tm = tm self.eigenvalues, self.eigenvectors = cm.diagonalize() # : array: transition quantities coming from the multipoles self.transitions = tm.get_transitions() scpr = np.array(np.dot(self.eigenvectors, self.transitions)) #: array: oscillator strenghts components of the transitions defined # as the square of $\int w_a(\mathbf r) r_i $ self.oscillator_strenghts = np.array([t**2 for t in scpr[:, 0:3]]) # : array: average of all the components of the OS self.avg_os = np.average(self.oscillator_strenghts, axis=1) self.alpha_prime = 2.0*self.oscillator_strenghts / \ self.eigenvalues[:, np.newaxis] """ array: elements of the integrand of the statical polarizability in the space of the excitations """ self._indices_for_spin_comparison = \ self._get_indices_for_spin_comparison() self.identify_singlet_and_triplets(1.e-5) def _get_indices_for_spin_comparison(self): inds = [[], []] inds0 = [] # get the indices for comparison, take the minimum between the spins if len(self.cm.norb_occ) == 1: nocc = self.cm.norb_occ[0] nvirt = self.cm.norb_virt[0] nos = [nocc, nocc] nvs = [nvirt, nvirt] else: nocc = min(self.cm.norb_occ) nvirt = min(self.cm.norb_virt) nos = self.cm.norb_occ nvs = self.cm.norb_virt for ispin in [0, 1]: for a in range(nvirt): for p in range(nocc): inds[ispin].append([p, a, ispin]) for a in range(nvirt, nvs[ispin]): for p in range(nocc, nos[ispin]): inds0.append([p, a, ispin]) transA = transition_indexes( self.cm.norb_occ, self.cm.norb_virt, inds[0]) transB = transition_indexes( self.cm.norb_occ, self.cm.norb_virt, inds[1]) trans0 = transition_indexes(self.cm.norb_occ, self.cm.norb_virt, inds0) return transA, transB, trans0
[docs] def spectrum_curves(self, omega, slice=None, weights=None): """Calculate spectrum curves. Provide the set of the curves associated to the weights. The resulting curves might then be used to draw the excitation spectra. Args: omega (array): the linspace used for the plotting, of shape ``(n,)``. Must be provided in Atomic Units slice (array): the lookup array that has to be considered. if Not provided the entire range is assumed weights (array): the set of arrays used to weight the spectra. Must have shape ``(rank,m)``, where ``rank`` is equal to the number of eigenvalues. If m is absent it is assumed to be 1. When not specified, it defaults to the average oscillator strenghts. Returns: array: a set of spectrum curves, of shape equal to ``(n,m)``, where ``n`` is the shape of ``omega`` and ``m`` the size of the second dimension of ``weights``. """ if slice is None: oo = self.eigenvalues[:, np.newaxis] - omega**2 wgts = weights if weights is not None else self.avg_os else: oo = self.eigenvalues[slice, np.newaxis] - omega**2 oo = oo[0] wgts = weights if weights is not None else self.avg_os[slice] return np.dot(2.0/oo.T, wgts)
[docs] def identify_singlet_and_triplets(self, tol=1.e-5): """ Find the lookup tables that select the singlets and the triplets among the excitations Args: tol (float): tolerance to be applied to recognise the spin character """ sings = [] trips = [] for exc in range(len(self.eigenvalues)): sing, trip = self.project_on_spin(exc, tol) if sing: sings.append(exc) if trip: trips.append(exc) if len(sings) > 0: self.singlets = (np.array(sings),) """array: lookup table of the singlet excitations""" if len(trips) > 0: self.triplets = (np.array(trips),) """array: lookup table of the triplet excitations"""
def _project_on_occ(self, exc): """ Project a given eigenvector on the occupied orbitals. In the spin averaged case consider all the spin indices nonetheless """ norb_occ = self.cm.norb_occ norb_virt = self.cm.norb_virt pProj_spin = [] for ispin, norb in enumerate(norb_occ): pProj = np.zeros(norb) for iorb in range(norb): harvest = _collection_indexes_iocc( iorb, self.cm.norb_virt, spin=None if len(norb_occ) == 1 else ispin) inds = np.array(transition_indexes( norb_occ, norb_virt, harvest)) pProj[iorb] = np.sum(np.ravel(self.eigenvectors[exc, inds])**2) pProj_spin.append(pProj) return pProj_spin
[docs] def project_on_spin(self, exc, tol=1.e-8): """ Control if an excitation has a Singlet or Triplet character Args: exc (int): index of the excitation to be controlled Returns: tuple (bool,bool): couple of booleans indicating if the excitation is a singlet or a triplet respectively """ A, B, zero = [np.ravel(self.eigenvectors[exc, ind]) for ind in self._indices_for_spin_comparison] issinglet = np.linalg.norm(A-B) < tol istriplet = np.linalg.norm(A+B) < tol return issinglet, istriplet
# print (self.eigenvalues[exc], np.linalg.norm(A), np.linalg.norm(B), # A-B,A+B, np.linalg.norm(zero)) def _get_threshold(self, pProj_spin, th_energies, tol): """ Identify the energy which is associated to the threshold of a given excitation. The tolerance is used to discriminate the component """ ths = -1.e100 for proj, en in zip(pProj_spin, th_energies): norb = len(en) pProj = proj.tolist() pProj.reverse() imax = norb-1 for val in pProj: if val > tol: break imax -= 1 ths = max(ths, en[imax]) return ths
[docs] def split_excitations(self, evals, tol, nexc='all'): """Separate the excitations in channels. This methods classify the excitations according to the channel they belong, and determine if a given excitation might be considered as a belonging to a discrete part of the spectrum or not. Args: evals (BandArray): the eigenvalues as they are provided (for instance) from a `Logfile` class instance. tol (float): tolerance for determining the threshold nexc (int,str): number of excitations to be analyzed. If ``'all'`` then the entire set of excitations are analyzed. """ self.determine_occ_energies(evals) self.identify_thresholds(self.occ_energies, tol, len( self.eigenvalues) if nexc == 'all' else nexc)
[docs] def identify_thresholds(self, occ_energies, tol, nexc): """Identify the thresholds per excitation. For each of the first ``nexc`` excitations, identify the energy value of its corresponding threshold. This value is determined by projecting the excitation components on the occupied states and verifying that their norm for the highest energy level is below a given tolerance. Args: occ_energies (tuple of array-like): contains the list of the energies of the occupied states per spin channel tol (float): tolerance for determining the threshold nexc (int): number of excitations to be analyzed """ # : Norm of the $w_p^a$ states associated to each excitation self.wp_norms = [] threshold_energies = [] for exc in range(nexc): proj = self._project_on_occ(exc) self.wp_norms.append(proj) threshold_energies.append( self._get_threshold(proj, occ_energies, tol)) # : list: identified threshold for inspected excitations self.threshold_energies = np.array(threshold_energies) self.excitations_below_threshold = np.where( np.abs(self.threshold_energies) >= np.sqrt( self.eigenvalues[0:len(self.threshold_energies)])) """ array: Indices of the excitations which lie below their corresponding threshold """ self.excitations_above_threshold = np.where( np.abs(self.threshold_energies) < np.sqrt(self.eigenvalues[0:len(self.threshold_energies)])) """ array: Indices of the excitations which lie above their corresponding threshold """
[docs] def determine_occ_energies(self, evals): """ Extract the occupied energy levels from a Logfile BandArray structure, provided the tuple of the number of occupied states Args: evals (BandArray): the eigenvalues as they are provided (for instance) from a `Logfile` class instance. """ norb_occ = self.cm.norb_occ occ_energies = [] # istart=0 for ispin, norb in enumerate(norb_occ): # range(len(norb_occ)): # istart:istart+norb_occ[ispin]])) occ_energies.append(np.array(evals[0][ispin][0:norb])) # istart+=norb_tot[ispin] # : array: energies of the occupied states out of the logfile self.occ_energies = occ_energies # : float: lowest threshold of the excitations. All excitations are # discrete below this level self.first_threshold = abs( max(np.max(self.occ_energies[0]), np.max(self.occ_energies[-1])))
[docs] def plot_alpha(self, **kwargs): """Plot the imaginary part. Plot the real or imaginary part of the dynamical polarizability. Keyword Arguments: real (bool): True if real part has to be plotted. The imaginary part is plotted otherwise eta (float): Value of the complex imaginary part. Defaults to 1.e-2. group (str): see :meth:`lookup` **kwargs: other arguments that might be passed to the :meth:`plot` method of the :mod:`matplotlib.pyplot` module. Returns: :mod:`matplotlib.pyplot`: the reference to :mod:`matplotlib.pyplot` module. """ import matplotlib.pyplot as plt from futile.Utils import kw_pop emax = np.max(np.sqrt(self.eigenvalues))*HaeV kwargs, real = kw_pop('real', False, **kwargs) plt.xlim(xmax=emax) if real: plt.ylabel(r'$\mathrm{Re} \alpha$ (AU)', size=14) else: plt.ylabel(r'$\mathrm{Im} \alpha$', size=14) plt.yticks([]) plt.xlabel(r'$\omega$ (eV)', size=14) if hasattr(self, 'first_threshold'): eps_h = self.first_threshold*HaeV plt.axvline(x=eps_h, color='black', linestyle='--') kwargs, eta = kw_pop('eta', 1.e-2, **kwargs) omega = np.linspace(0.0, emax, 5000)+2.0*eta*1j kwargs, group = kw_pop('group', 'all', **kwargs) slice = self.lookup(group) spectrum = self.spectrum_curves(omega, slice=slice) toplt = spectrum.real if real else spectrum.imag pltkwargs = dict(c='black', linewidth=1.5) pltkwargs.update(kwargs) plt.plot(omega*HaeV, toplt, **pltkwargs) return plt
[docs] def lookup(self, group): """ Identify the group of the excitations according to the argument Args: group (str): A string chosen between * ``"all"`` : provides the entire set of excitations (:py:class:`None` instead of the lookup array) * ``"bt"`` : provides only the excitations below threshold * ``"at"`` : provides only the excitations above threshold * ``"singlets"`` : provides the index of the excitations that have a singlet character * ``"triplets"`` : provides the index of the excitations that have a triplet character """ slice = None if group == 'bt': slice = self.excitations_below_threshold if group == 'at': slice = self.excitations_above_threshold if group == 'singlets': slice = self.singlets if group == 'triplets': slice = self.triplets return slice
[docs] def plot_excitation_landscape(self, **kwargs): """ Represent the excitation landscape as splitted in the excitations class Args: **kwargs: keyword arguments to be passed to the `pyplot` instance. The ``xlabel``, ``ylabel`` as well as ``xlim`` are already set. Returns: :mod:`matplotlib.pyplot`: the reference to :mod:`matplotlib.pyplot` module. Example: >>> ex=Excitations(cm,tm) >>> ex.split_excitations(evals=...,tol=1.e-4,nexc=...) >>> ex.plot_excitation_landscape(title='Excitation landscape') """ import matplotlib.pyplot as plt Emin = 0.0 Emax = np.max(np.sqrt(self.eigenvalues))*HaeV for level in self.occ_energies[0]: eng_th = level*HaeV plt.plot((Emin, eng_th), (level, level), '--', c='red', linewidth=1) plt.plot((eng_th, Emax), (level, level), '-', c='red', linewidth=1) plt.scatter(abs(eng_th), level, marker='x', c='red') ind_bt = self.excitations_below_threshold exc_bt = np.sqrt(self.eigenvalues)[ind_bt] lev_bt = self.threshold_energies[ind_bt] plt.scatter(HaeV*exc_bt, lev_bt, s=16, marker='o', c='black') ind_at = self.excitations_above_threshold exc_at = np.sqrt(self.eigenvalues)[ind_at] lev_at = self.threshold_energies[ind_at] plt.scatter(HaeV*exc_at, lev_at, s=14, marker='s', c='blue') plt.xlabel('energy (eV)') plt.ylabel('Threshold energy (Ha)') plt.xlim(xmin=Emin-1, xmax=Emax) for attr, val in kwargs.items(): if type(val) == dict: getattr(plt, attr)(**val) else: getattr(plt, attr)(val) return plt
[docs] def dos_dict(self, group='all'): """Dictionary for DoS creation. Creates the keyword arguments that have to be passed to the `meth:BigDFT.DoS.append` method of the `DoS` class Args: group (str): see :meth:`lookup` Returns: :py:class:`dict`: kewyord arguments that can be passed to the `meth:BigDFT.DoS.append` method of the :class:`DoS.DoS` class """ ev = np.sqrt(self.eigenvalues) slice = self.lookup(group) if slice is not None: ev = ev[slice] return dict(energies=np.array([np.ravel(ev)]), units='AU')
[docs] def dos(self, group='all', **kwargs): """Density of States of the Excitations. Provides an instance of the :class:`~BigDFT.DoS.DoS` class, corresponding to the Excitations instance. Args: group (str): see :meth:`lookup` **kwargs: other arguments that might be passed to the :class:`DoS.DoS` instantiation Returns: :class:`DoS.DoS`: instance of the Density of States class """ from BigDFT.DoS import DoS kwa = self.dos_dict(group=group) kwa['energies'] = kwa['energies'][0] if hasattr(self, 'first_threshold'): kwa['fermi_level'] = self.first_threshold else: kwa['fermi_level'] = 0.0 kwa.update(kwargs) return DoS(**kwa)
[docs] def plot_Sminustwo(self, coord, alpha_ref=None, group='all'): """Inspect S-2 sum rule. Provides an handle to the plotting of the $S^{-2}$ sum rule, which should provide reference values for the static polarizability tensor. Args: coord (str): the coordinate used for inspection. May be ``'x'``, ``'y'`` or ``'z'``. alpha_ref (list): diagonal of the reference static polarizability tensor (for instance calculated via finite differences). If present the repartition of the contribution of the various groups of excitations is plotted. group (str): see :meth:`lookup` Returns: reference to :mod:`matplotlib.pyplot` module. """ import matplotlib.pyplot as plt idir = ['x', 'y', 'z'].index(coord) fig, ax1 = plt.subplots() ax1.set_xlabel('energy (eV)', size=14) plt.ylabel(r'$\alpha_{'+coord+coord+r'}$ (AU)', size=14) if alpha_ref is not None: plt.axhline(y=alpha_ref[idir], color='r', linestyle='--') if hasattr(self, 'first_threshold'): eps_h = abs(HaeV*self.first_threshold) plt.axvline(x=eps_h, color='black', linestyle='--') e = np.sqrt(self.eigenvalues)*HaeV w_ii = self.alpha_prime[:, idir] slice = self.lookup(group) if slice is not None: e = e[slice] w_ii = w_ii[slice] ax1.plot(e, np.cumsum(w_ii)) ax2 = ax1.twinx() ax2.plot(e, w_ii, color='grey', linestyle='-') plt.ylabel(r'$w_{'+coord+coord+r'}$ (AU)', size=14) return plt
[docs]def get_alpha_energy(log, norb, nalpha): return log.evals[0][0][norb+nalpha-1]
[docs]def identify_contributions(numOrb, na, exc, C_E2): pProj = np.zeros(numOrb*2) for p in range(numOrb): for spin in [0, 1]: # sum over all the virtual orbital and spin for alpha in range(na): # extract the value of the index of C_E2 elements = transition_indexes( [numOrb], [na], [[p, alpha, spin]]) for el in elements: pProj[p+numOrb*spin] += C_E2[exc][el]**2 pProj = pProj[0:numOrb]+pProj[numOrb:2*numOrb] # halves the components return pProj
[docs]def get_p_energy(log, norb): return log.evals[0][0][0:norb]
[docs]def get_threshold(pProj, th_energies, th_levels, tol): norb = len(th_energies) pProj = pProj.tolist() pProj.reverse() imax = norb-1 for val in pProj: if val > tol: break imax -= 1 return [th_levels[imax], th_energies[imax]]